"""
Phase 3: Level 2 Commodity-Demographics Deep Dive

Key questions:
1. WHY does Z₁ amplify in commodity economies? (savings reinforcement)
2. Does the development threshold differ for commodity vs non-commodity?
3. Does the doom loop operate differently with commodity revenue?
4. Is there a "double cliff" — commodity depletion + aging?
5. Gravity model: bilateral flows for commodity exporters
6. Does commodity dependence change the fiscal transmission of aging?
"""

import sys
import os
import pandas as pd
import numpy as np
from scipy import stats

sys.path.insert(0, '/mnt/c/demographics_capital_flows/multilateral/src')
from model import PanelGLS

OUT = '/mnt/c/demographics_capital_flows/eu_demographics/output/tables'
os.makedirs(OUT, exist_ok=True)

# ── Load data ──────────────────────────────────────────────────────────
print("Loading data...")
panel = pd.read_csv('/mnt/c/demographics_capital_flows/multilateral/140_country/data/processed/full_panel_with_resources.csv')
panel = panel[panel['year'] <= 2024].copy()

# Also load fiscal panel for doom loop variables
fiscal = pd.read_csv('/mnt/c/demographics_capital_flows/fiscal_dominance/data/processed/fiscal_panel.csv')

# Merge fiscal stress if available
if 'fiscal_stress' in fiscal.columns or 'fd_binary' in fiscal.columns:
    fisc_regimes = pd.read_csv('/mnt/c/demographics_capital_flows/fiscal_dominance/data/processed/fiscal_regimes.csv')
    panel = panel.merge(fisc_regimes[['iso3', 'year', 'fd_binary', 'fiscal_stress']],
                       on=['iso3', 'year'], how='left')

# Ensure resource variables exist
if 'resource_rents_gdp' not in panel.columns:
    print("ERROR: Run phase2 first.")
    sys.exit(1)

# Create additional variables
panel['commodity_exporter'] = (panel['resource_rents_gdp'] >= 10).astype(float)
panel['high_resource'] = (panel['resource_rents_gdp'] >= 5).astype(float)
panel['Z1_x_resource'] = panel['Z_1'] * panel['resource_rents_gdp']
panel['Z1_x_commodity'] = panel['Z_1'] * panel['commodity_exporter']

# Income interactions
if 'gdp_pc_ppp' not in panel.columns:
    for col in ['rgdpna', 'gdp_usd']:
        if col in panel.columns:
            panel['gdp_proxy'] = panel[col]
            break
else:
    panel['gdp_proxy'] = panel['gdp_pc_ppp']

controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth', 'log_rel_opw', 'kaopen']
gls = PanelGLS()


def run_model(name, dep, indep, data, quiet=False):
    """Run PanelGLS, print, return dict."""
    est = data[['iso3', 'year', dep] + indep].dropna()
    if len(est) < 50:
        print(f"  {name}: insufficient obs ({len(est)})")
        return None
    y = est[dep].values
    X = est[indep].values
    m = PanelGLS()
    m.fit(y, X, est['iso3'].values, est['year'].values)
    if not quiet:
        print(f"\n  {name}: N={m.n_obs}, R²={m.r_squared:.4f}")
        for i, v in enumerate(indep):
            s = '***' if m.pvalues[i]<0.01 else '**' if m.pvalues[i]<0.05 else '*' if m.pvalues[i]<0.1 else ''
            print(f"    {v:35s} {m.beta[i]:10.4f} ({m.se[i]:.4f}){s}")
    return {'model': name, 'n_obs': m.n_obs, 'r_squared': m.r_squared,
            'vars': indep, 'beta': m.beta, 'se': m.se, 'pvalues': m.pvalues}


# ═══════════════════════════════════════════════════════════════════════
# SECTION 1: WHY DOES Z₁ AMPLIFY IN COMMODITY ECONOMIES?
# Hypothesis: aging commodity exporters save windfalls rather than
# spending them → lifecycle savings + commodity savings reinforce
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 1: Savings Reinforcement Hypothesis")
print("=" * 70)

# Test: Z₁ → savings rate, moderated by commodity status
if 'gross_savings_gdp' in panel.columns:
    print("\n--- Savings channel ---")
    run_model('Savings: Baseline', 'gross_savings_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls, panel)

    run_model('Savings: + Resource Rents', 'gross_savings_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls + ['resource_rents_gdp'], panel)

    run_model('Savings: + Z₁ × Resource', 'gross_savings_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls + ['resource_rents_gdp', 'Z1_x_resource'], panel)

    # Split sample
    low_res = panel[panel['resource_rents_gdp'] < 5]
    high_res = panel[panel['resource_rents_gdp'] >= 5]
    comm = panel[panel['resource_rents_gdp'] >= 10]

    run_model('Savings: Low Resource (<5%)', 'gross_savings_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls, low_res)

    run_model('Savings: Commodity (≥10%)', 'gross_savings_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls, comm)

# Test: Z₁ → investment rate, moderated by commodity status
if 'gross_investment_gdp' in panel.columns:
    print("\n--- Investment channel ---")
    run_model('Investment: Baseline', 'gross_investment_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls, panel)

    run_model('Investment: + Z₁ × Resource', 'gross_investment_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls + ['resource_rents_gdp', 'Z1_x_resource'], panel)

    run_model('Investment: Low Resource', 'gross_investment_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls, low_res)

    run_model('Investment: Commodity (≥10%)', 'gross_investment_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls, comm)

# Savings-investment gap (= CA, approximately)
if 'savings_investment_gap' in panel.columns:
    print("\n--- S-I gap channel ---")
    run_model('S-I Gap: + Z₁ × Resource', 'savings_investment_gap',
              ['Z_1', 'Z_2', 'Z_3'] + controls + ['resource_rents_gdp', 'Z1_x_resource'], panel)

# ═══════════════════════════════════════════════════════════════════════
# SECTION 2: DEVELOPMENT THRESHOLD × COMMODITY
# Does income interact differently with demographics for commodity
# vs non-commodity economies?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 2: Development Threshold × Commodity Dependence")
print("=" * 70)

if 'gdp_proxy' in panel.columns and panel['gdp_proxy'].notna().any():
    # Create income bins
    panel['log_gdp'] = np.log(panel['gdp_proxy'].clip(lower=1))
    panel['Z1_x_log_gdp'] = panel['Z_1'] * panel['log_gdp']
    panel['Z1_x_log_gdp_x_comm'] = panel['Z_1'] * panel['log_gdp'] * panel['commodity_exporter']
    panel['log_gdp_x_comm'] = panel['log_gdp'] * panel['commodity_exporter']

    print("\n--- Triple interaction: Z₁ × Income × Commodity ---")
    run_model('Threshold: Z₁ × log(GDP)', 'ca_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls + ['log_gdp', 'Z1_x_log_gdp'], panel)

    run_model('Threshold: + Commodity triple', 'ca_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + controls +
              ['log_gdp', 'Z1_x_log_gdp', 'commodity_exporter',
               'Z1_x_commodity', 'log_gdp_x_comm', 'Z1_x_log_gdp_x_comm'], panel)

    # Split by income tier and commodity status
    if panel['gdp_proxy'].notna().sum() > 100:
        med_gdp = panel['gdp_proxy'].median()
        low_inc_noncomm = panel[(panel['gdp_proxy'] < med_gdp) & (panel['commodity_exporter'] == 0)]
        low_inc_comm = panel[(panel['gdp_proxy'] < med_gdp) & (panel['commodity_exporter'] == 1)]
        high_inc_noncomm = panel[(panel['gdp_proxy'] >= med_gdp) & (panel['commodity_exporter'] == 0)]
        high_inc_comm = panel[(panel['gdp_proxy'] >= med_gdp) & (panel['commodity_exporter'] == 1)]

        print("\n--- Split by income × commodity status ---")
        for label, sub in [('Low-inc Non-commodity', low_inc_noncomm),
                           ('Low-inc Commodity', low_inc_comm),
                           ('High-inc Non-commodity', high_inc_noncomm),
                           ('High-inc Commodity', high_inc_comm)]:
            run_model(f'Threshold: {label}', 'ca_gdp',
                     ['Z_1', 'Z_2', 'Z_3'] + controls, sub)

# ═══════════════════════════════════════════════════════════════════════
# SECTION 3: FISCAL TRANSMISSION — DOES COMMODITY REVENUE CHANGE
# HOW DEMOGRAPHICS AFFECT FISCAL BALANCES?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 3: Fiscal Transmission × Commodity")
print("=" * 70)

# Demographics → fiscal balance, interacted with commodity
tri_controls = ['nfa_gdp_lag', 'rgdp_growth']  # avoid fiscal_bal on RHS

run_model('Fiscal: Z₁ → fiscal_bal (baseline)', 'fiscal_bal_gdp',
          ['Z_1', 'Z_2', 'Z_3'] + tri_controls, panel)

run_model('Fiscal: + Z₁ × Resource', 'fiscal_bal_gdp',
          ['Z_1', 'Z_2', 'Z_3'] + tri_controls + ['resource_rents_gdp', 'Z1_x_resource'], panel)

# Split sample
low_res = panel[panel['resource_rents_gdp'] < 5]
high_res = panel[panel['resource_rents_gdp'] >= 5]
comm = panel[panel['resource_rents_gdp'] >= 10]

run_model('Fiscal: Low Resource (<5%)', 'fiscal_bal_gdp',
          ['Z_1', 'Z_2', 'Z_3'] + tri_controls, low_res)

run_model('Fiscal: Commodity (≥10%)', 'fiscal_bal_gdp',
          ['Z_1', 'Z_2', 'Z_3'] + tri_controls, comm)

# Demographics → government debt, interacted
if 'govt_debt_gdp' in panel.columns:
    print("\n--- Debt channel ---")
    run_model('Debt: Z₁ → debt/GDP', 'govt_debt_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + tri_controls, panel)

    run_model('Debt: + Z₁ × Resource', 'govt_debt_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + tri_controls + ['resource_rents_gdp', 'Z1_x_resource'], panel)

    run_model('Debt: Commodity (≥10%)', 'govt_debt_gdp',
              ['Z_1', 'Z_2', 'Z_3'] + tri_controls, comm)

# ═══════════════════════════════════════════════════════════════════════
# SECTION 4: SOVEREIGN SPREADS × COMMODITY
# Do markets price demographic risk differently for commodity exporters?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 4: Sovereign Spreads × Commodity")
print("=" * 70)

spread = pd.read_csv('/mnt/c/demographics_capital_flows/sovereign_spreads/data/processed/spread_panel.csv')
# Merge resource rents
rents = pd.read_csv('/mnt/c/demographics_capital_flows/multilateral/140_country/data/raw/wdi_resource_rents.csv')
spread = spread.merge(rents[['iso3', 'year', 'resource_rents_gdp']], on=['iso3', 'year'], how='left')
spread['commodity_exporter'] = (spread['resource_rents_gdp'] >= 10).astype(float)
spread['Z1_x_resource'] = spread['Z_1'] * spread['resource_rents_gdp']
spread['Z1_x_commodity'] = spread['Z_1'] * spread['commodity_exporter']

spread_dep = None
for col in ['sovereign_spread', 'spread_vs_us', 'govt_bond_10y']:
    if col in spread.columns and spread[col].notna().sum() > 200:
        spread_dep = col
        break

if spread_dep:
    spread_controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']
    spread_est = spread[spread['year'] <= 2024]

    run_model(f'Spread: Z₁ → {spread_dep}', spread_dep,
              ['Z_1', 'Z_2', 'Z_3'] + spread_controls, spread_est)

    run_model(f'Spread: + Z₁ × Resource', spread_dep,
              ['Z_1', 'Z_2', 'Z_3'] + spread_controls + ['resource_rents_gdp', 'Z1_x_resource'],
              spread_est)

    run_model(f'Spread: + Z₁ × Commodity', spread_dep,
              ['Z_1', 'Z_2', 'Z_3'] + spread_controls + ['commodity_exporter', 'Z1_x_commodity'],
              spread_est)

    # Split
    spread_noncomm = spread_est[spread_est['commodity_exporter'] == 0]
    spread_comm = spread_est[spread_est['commodity_exporter'] == 1]
    run_model(f'Spread: Non-Commodity', spread_dep,
              ['Z_1', 'Z_2', 'Z_3'] + spread_controls, spread_noncomm)
    run_model(f'Spread: Commodity', spread_dep,
              ['Z_1', 'Z_2', 'Z_3'] + spread_controls, spread_comm)

# ═══════════════════════════════════════════════════════════════════════
# SECTION 5: THE DOUBLE CLIFF — AGING + COMMODITY DEPLETION
# Which commodity exporters are aging fast enough that they face
# both demographic pressure AND eventual commodity revenue decline?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 5: The Double Cliff — Aging × Commodity Depletion")
print("=" * 70)

# Identify commodity exporters with aging trajectories
# Use current resource rents + demographic projections
recent = panel[(panel['year'] >= 2018) & (panel['year'] <= 2021)]
country_snap = recent.groupby('iso3').agg(
    mean_rents=('resource_rents_gdp', 'mean'),
    Z1_current=('Z_1', 'last'),
    oadr_current=('old_dep', 'last'),
    ca_current=('ca_gdp', 'last'),
    nfa_current=('nfa_gdp', 'last')
).dropna(subset=['mean_rents'])

# Get projections
full_proj = pd.read_csv('/mnt/c/demographics_capital_flows/multilateral/140_country/data/processed/full_panel.csv')
proj_2050 = full_proj[full_proj['year'] == 2050][['iso3', 'Z_1', 'old_dep']].rename(
    columns={'Z_1': 'Z1_2050', 'old_dep': 'oadr_2050'})

country_snap = country_snap.merge(proj_2050, left_index=True, right_on='iso3')
country_snap['Z1_change'] = country_snap['Z1_2050'] - country_snap['Z1_current']
country_snap['oadr_change'] = country_snap['oadr_2050'] - country_snap['oadr_current']

# Filter to commodity exporters (≥10% rents)
comm_aging = country_snap[country_snap['mean_rents'] >= 10].sort_values('Z1_change', ascending=False)

print("\nCommodity exporters (≥10% rents) ranked by aging speed (Z₁ change to 2050):")
print(f"{'Country':6s} {'Rents%':>7s} {'Z₁ now':>7s} {'Z₁ 2050':>8s} {'ΔZ₁':>6s} {'OADR now':>9s} {'OADR 2050':>10s} {'CA':>6s}")
print("-" * 70)
for _, row in comm_aging.iterrows():
    print(f"{row['iso3']:6s} {row['mean_rents']:7.1f} {row['Z1_current']:7.3f} "
          f"{row['Z1_2050']:8.3f} {row['Z1_change']:6.3f} "
          f"{row['oadr_current']*100:8.1f}% {row['oadr_2050']*100:9.1f}% "
          f"{row['ca_current']:6.1f}")

# Identify "double cliff" countries: high rents AND fast aging
double_cliff = comm_aging[comm_aging['Z1_change'] > 0.5]
print(f"\n'Double cliff' countries (commodity exporter + Z₁ rising >0.5 by 2050):")
print(f"  {', '.join(double_cliff['iso3'].tolist())}")

# ═══════════════════════════════════════════════════════════════════════
# SECTION 6: COMMODITY ECONOMIES AND CRISIS RISK
# Does resource dependence change the demographic EWS signal?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 6: Crisis Risk × Commodity")
print("=" * 70)

# Check for crisis variables
crisis_cols = [c for c in panel.columns if 'crisis' in c.lower() or 'banking' in c.lower()]
if crisis_cols:
    print(f"  Crisis variables found: {crisis_cols}")
    for cc in crisis_cols[:2]:
        if panel[cc].notna().sum() > 200:
            run_model(f'Crisis: Z₁ → {cc}', cc,
                     ['Z_1', 'Z_2', 'Z_3'] + controls, panel)
            run_model(f'Crisis: + Z₁ × Resource', cc,
                     ['Z_1', 'Z_2', 'Z_3'] + controls + ['resource_rents_gdp', 'Z1_x_resource'],
                     panel)
else:
    print("  No crisis variables in panel. Checking for fiscal stress...")

if 'fd_binary' in panel.columns:
    run_model('Fiscal Dominance: Z₁ → fd_binary', 'fd_binary',
              ['Z_1', 'Z_2', 'Z_3', 'nfa_gdp_lag', 'rgdp_growth'], panel)
    run_model('Fiscal Dominance: + Z₁ × Resource', 'fd_binary',
              ['Z_1', 'Z_2', 'Z_3', 'nfa_gdp_lag', 'rgdp_growth',
               'resource_rents_gdp', 'Z1_x_resource'], panel)

# ═══════════════════════════════════════════════════════════════════════
# SECTION 7: BILATERAL FLOWS — DO COMMODITY EXPORTERS ATTRACT
# DIFFERENT CAPITAL?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 7: Bilateral Flows to/from Commodity Exporters")
print("=" * 70)

bilateral = pd.read_csv('/mnt/c/demographics_capital_flows/gravity_bilateral/data/processed/bilateral_panel.csv')

# Add resource rents for destination
bilateral = bilateral.merge(
    rents[['iso3', 'year', 'resource_rents_gdp']].rename(
        columns={'iso3': 'iso_d', 'resource_rents_gdp': 'resource_rents_d'}),
    on=['iso_d', 'year'], how='left')

# Add for origin
bilateral = bilateral.merge(
    rents[['iso3', 'year', 'resource_rents_gdp']].rename(
        columns={'iso3': 'iso_o', 'resource_rents_gdp': 'resource_rents_o'}),
    on=['iso_o', 'year'], how='left')

bilateral['comm_dest'] = (bilateral['resource_rents_d'] >= 10).astype(float)
bilateral['comm_origin'] = (bilateral['resource_rents_o'] >= 10).astype(float)

if 'dZ_1' in bilateral.columns:
    bilateral['dZ1_x_comm_dest'] = bilateral['dZ_1'] * bilateral['comm_dest']
    bilateral['dZ1_x_comm_origin'] = bilateral['dZ_1'] * bilateral['comm_origin']

flow_vars = ['portfolio_total', 'portfolio_debt', 'portfolio_equity', 'fdi_outward']
gravity_controls = ['log_dist', 'contiguity', 'common_lang_official']
grav_avail = [c for c in gravity_controls if c in bilateral.columns]

for fv in flow_vars:
    if fv in bilateral.columns and bilateral[fv].notna().sum() > 500:
        bilateral[f'log_{fv}'] = np.log(bilateral[fv].clip(lower=1))

        print(f"\n--- {fv} ---")
        run_model(f'Gravity {fv}: dZ₁ baseline', f'log_{fv}',
                 ['dZ_1', 'dZ_2', 'dZ_3'] + grav_avail, bilateral)

        run_model(f'Gravity {fv}: + dZ₁ × comm_dest', f'log_{fv}',
                 ['dZ_1', 'dZ_2', 'dZ_3'] + grav_avail +
                 ['comm_dest', 'dZ1_x_comm_dest'], bilateral)

        run_model(f'Gravity {fv}: + dZ₁ × comm_origin', f'log_{fv}',
                 ['dZ_1', 'dZ_2', 'dZ_3'] + grav_avail +
                 ['comm_origin', 'dZ1_x_comm_origin'], bilateral)

        # Only run first two flow types to keep output manageable
        if fv == 'portfolio_equity':
            break

# ═══════════════════════════════════════════════════════════════════════
# SECTION 8: TIME VARIATION — HAS THE COMMODITY-DEMOGRAPHIC
# INTERACTION CHANGED OVER TIME?
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 8: Time Variation in Commodity-Demographic Interaction")
print("=" * 70)

# Roll through decades
base_vars = ['Z_1', 'Z_2', 'Z_3'] + controls
for period, (y1, y2) in [('1980-1999', (1980, 1999)),
                          ('2000-2009', (2000, 2009)),
                          ('2010-2021', (2010, 2021))]:
    sub = panel[(panel['year'] >= y1) & (panel['year'] <= y2)]
    r = run_model(f'Period {period}: + Z₁ × Resource', 'ca_gdp',
                  base_vars + ['resource_rents_gdp', 'Z1_x_resource'], sub, quiet=True)
    if r:
        z1_idx = r['vars'].index('Z_1')
        int_idx = r['vars'].index('Z1_x_resource')
        rr_idx = r['vars'].index('resource_rents_gdp')
        z1_s = '***' if r['pvalues'][z1_idx]<0.01 else '**' if r['pvalues'][z1_idx]<0.05 else '*' if r['pvalues'][z1_idx]<0.1 else ''
        int_s = '***' if r['pvalues'][int_idx]<0.01 else '**' if r['pvalues'][int_idx]<0.05 else '*' if r['pvalues'][int_idx]<0.1 else ''
        rr_s = '***' if r['pvalues'][rr_idx]<0.01 else '**' if r['pvalues'][rr_idx]<0.05 else '*' if r['pvalues'][rr_idx]<0.1 else ''
        print(f"  {period}: N={r['n_obs']}, R²={r['r_squared']:.4f}, "
              f"Z₁={r['beta'][z1_idx]:.3f}{z1_s}, "
              f"Z₁×Res={r['beta'][int_idx]:.4f}{int_s}, "
              f"Rents={r['beta'][rr_idx]:.3f}{rr_s}")

# ═══════════════════════════════════════════════════════════════════════
# SECTION 9: COMPREHENSIVE COUNTRY CLASSIFICATION
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("SECTION 9: Country Classification — Demographic-Commodity Matrix")
print("=" * 70)

# 2x2: Old/Young × Commodity/Non-commodity
recent_panel = panel[(panel['year'] >= 2015) & (panel['year'] <= 2021)]
country_class = recent_panel.groupby('iso3').agg(
    mean_Z1=('Z_1', 'mean'),
    mean_rents=('resource_rents_gdp', 'mean'),
    mean_ca=('ca_gdp', 'mean'),
    mean_nfa=('nfa_gdp', 'mean'),
    mean_fiscal=('fiscal_bal_gdp', 'mean'),
    mean_savings=('gross_savings_gdp', 'mean') if 'gross_savings_gdp' in recent_panel.columns else ('ca_gdp', 'count')
).dropna(subset=['mean_Z1', 'mean_rents'])

# Classify
country_class['old'] = country_class['mean_Z1'] > 0
country_class['commodity'] = country_class['mean_rents'] >= 10

quadrant_map = {
    (True, True): 'Old Commodity (Gulf/Russia)',
    (True, False): 'Old Non-Commodity (OECD core)',
    (False, True): 'Young Commodity (SSA oil/minerals)',
    (False, False): 'Young Non-Commodity (EM diversified)'
}

country_class['quadrant'] = country_class.apply(
    lambda r: quadrant_map[(r['old'], r['commodity'])], axis=1)

print("\n2×2 Classification:")
for q in quadrant_map.values():
    sub = country_class[country_class['quadrant'] == q]
    print(f"\n  {q}: {len(sub)} countries")
    print(f"    Mean Z₁: {sub['mean_Z1'].mean():.3f}, Mean CA: {sub['mean_ca'].mean():.2f}, "
          f"Mean Rents: {sub['mean_rents'].mean():.1f}%")
    if len(sub) <= 15:
        print(f"    Countries: {', '.join(sub.index.tolist())}")
    else:
        top5 = sub.nlargest(5, 'mean_rents')
        print(f"    Top 5 by rents: {', '.join(top5.index.tolist())}")

# Run CA regression for each quadrant
print("\n--- Quadrant-specific Z₁ effects on CA ---")
for q in quadrant_map.values():
    countries = country_class[country_class['quadrant'] == q].index.tolist()
    sub = panel[panel['iso3'].isin(countries)]
    r = run_model(f'  {q[:25]}', 'ca_gdp', base_vars, sub, quiet=True)
    if r:
        z1_idx = r['vars'].index('Z_1')
        s = '***' if r['pvalues'][z1_idx]<0.01 else '**' if r['pvalues'][z1_idx]<0.05 else '*' if r['pvalues'][z1_idx]<0.1 else ''
        print(f"  {q:40s}: N={r['n_obs']:5d}, Z₁={r['beta'][z1_idx]:8.3f}{s:3s}, R²={r['r_squared']:.4f}")

# ═══════════════════════════════════════════════════════════════════════
# SECTION 10: WRITE COMPREHENSIVE OUTPUT TABLE
# ═══════════════════════════════════════════════════════════════════════
print("\n" + "=" * 70)
print("Writing output tables...")
print("=" * 70)

# Save classification
country_class.to_csv(os.path.join(OUT, 'commodity_demographic_classification.csv'))

# Save double cliff
comm_aging.to_csv(os.path.join(OUT, 'double_cliff_risk.csv'), index=False)

print("  Written: commodity_demographic_classification.csv")
print("  Written: double_cliff_risk.csv")
print("\n✓ Phase 3 complete.")
